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Abstract We propose a novel optimization algorithm for continuous functions using 
geodesics and contours under conformal mapping. The algorithm can find multiple 
optima by first following a geodesic curve to a local optimum then traveling to the 
next search area by following a contour curve. To improve the efficiency, Newton- 
Raphson algorithm is also employed in local search steps. A proposed jumping mech¬ 
anism based on realized geodesics enables the algorithm to jump to a nearby region 
and consequently avoid trapping in local optima. Conformal mapping is used to re¬ 
solve numerical instability associated with solving the classical geodesic equations. 
Geodesic flows under conformal mapping are constructed numerically by using lo¬ 
cal quadratic approximation. The parameters in the algorithm are adaptively cho¬ 
sen to reflect local geometric features of the objective function. Comparisons with 
many commonly used optimization algorithms including gradient, trust region, ge¬ 
netic algorithm and global search methods have shown that the proposed algorithm 
outperforms most widely used methods in almost all test cases with only a couple of 
exceptions. 
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1 Introduction 

Optimization is an essential process in scientific investigation. There are many effec¬ 
tive and efficient methods proposed in the literature, see for example (D, 12, II and 
g). We propose a new algorithm that solves optimization problems from a point of 
view that is different than most of the major methods in the literature. The proposed 
method builds one dimensional paths or curves and travels on it with a constant speed 
to search for the global optimum. 

The main idea of the proposed optimization algorithm is built upon geodesics 
which are a generalization of straight lines in Euclidean space that minimizes the 
non-Euclidean distance between two points on a given manifold defined by the ob¬ 
jective function. The paths of optimization are constructed numerically by using a 
local quadratic approximation since there is no analytical solution to the geodesic 
equation in a general manifold. To avoid numerical instabilities in converting ill- 
conditioned matrices, the geodesics are constructed on a manifold under conformal 
mapping which preserves intrinsic geometrical features of the objective function. The 
algorithm will then follow either the geodesics or contours under conformal mapping 
to search for the optimum. Each contour curve could provide a bridge to a new search 
area and the constant speed enforced by the algorithm ensures that the search never 
stops or be trapped at a stationary point. Even with carefully constructed geodesics, 
the proposal algorithm can still be trapped within one region. In order to search an¬ 
other promising region, we also build a jumping mechanism by examining the values 
of the objective function along the geodesics to detect any potential and hidden influ¬ 
ence to the geodesic flow from an nearby optimum. 

The algorithm can be further improved by integrating with other search meth¬ 
ods. For example, one can use a few points along the geodesic as starting points 
to a Quasi-Newton algorithm to improve computational efficiency. From the Quasi- 
Newton outputs the algorithm is able to change its parameters adaptively for oscil¬ 
lating and smooth objective functions. We found that the resulting algorithm per¬ 
forms well in both types of functions in moderately high dimensions. Furthermore, 
we built a stopping criterion for the algorithm using Quasi-Newton methods by set¬ 
ting a threshold on the number of the maximum found within tolerance. 

The remaining part of the paper is organized as follows. In Section[2]we give an 
introductory review on differential geometry. Theoretical properties of the proposed 
algorithm are established in Section [3] In Section [TJ we give a general description 
of the algorithm and we also describe the method of choosing the parameters adap¬ 
tively. Numerical results comparing the proposed algorithm with the Quasi-Newton, 
genetic algorithm, wedge trust region methods and the global search function in Mat- 
lab’s global optimization package are provided in Section [5] Finally, the conclusion 
is given in Section[6] 


2 The Main Idea 

We generalize the line search method with geodesics in order to discover multiple 
maxima on a manifold conformally related to R". 
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2.1 Geodesics and Geodesic Equations 


We consider a topological manifold that is a second countable and locally compact 
Hausdorff space. It is also connected and completely regular. Detailed discussions 
can be found in Q and 0. A Riemannian metric on a smooth and differentiable 
manifold M is a 2-tensor field T 2 (M) that is symmetric and positive definite. A 
Riemannian metric thus determines an inner product on each tangent space T p {M), 
which is typically written as g(U, V") for U, V £ T p (M). For an Euclidean space, the 
metric matrix (or just metric henceforth) in component form is the Kronecka delta, 
i.e gij = Sij. The inner product g(U, V) in Euclidean reduces to the dot product, 
JA • SijCPVJ, where the sum is over all dimensions. In the Einstein summation con¬ 
vention, it is understood that repeated indices are summed over and the inner product 
is expressed as g l jU l V 3 . 

A geodesic is defined to be the path of minimum length for two given distinct 
points in a connected manifold. It is simply a straight line in Euclidean space. In 
a non-flat manifold, however, it is a curve and no longer a straight line. Let X l (t) 
denote the local coordinate for the i-th dimension for a parameter t which is a time 
step in our case. The geodesic is then characterized by a set of partial differential 
equations, using the Einstein summation convention: 


d 2 X'{t) , dX*(t) dX k {t) 

dt 2 dt dt 


(1) 


where I '- k are Christoffel symbols defined to be 
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g^ is the metric and g’ :l is the inverse metric. 

There exists a unique vector field on the tangent bundle of manifold, denoted 
as TM, whose trajectories are of the form (7(f), 7 (i)) where 7 is the geodesic. 
Geodesics play an important role in General Relativity, see 0, where the path of a 
planet orbiting around a star is the projection of a geodesic of the curved 4-D space- 
time geometry around the star onto a 3-D Euclidean space. 


2.2 Conformal Mapping 

Numerical calculations of the Christoffel symbols involve the inversion of the metric 
and can be unstable and computationally costly. One strategy to avoid such com¬ 
plications is to calculate the Christoffel symbols in a manifold where the metric is 
easily inverted and then map the results to the manifold desired. We consider the case 
where the manifold containing information about the objective function is mapped 
from R", where the metric and the inverse metric is the Kronecka delta. In such a 
case an analytic expression for the Christoffel symbols is available and the costly 
matrix inversion is avoided. 
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Each local neighborhood in the new manifold is holomorphic to R" . The re¬ 
sulting metric under the conformal mapping is said to be conformally related to the 
Euclidean metric 

9ij = 'P(x) 2 6 ij , 

where the scale factor &(x) = and <p(x) is a real valued objective function. 

Trivially, the manifolds obtained this way are Riemannian as the metric tensors are 
both symmetric and positive definite. The existence of such mappings is trivial since 
we assume the existence of cf>(x) to begin with. 


3 Theoretical Properties of the Geodesic in a Conformally Flat Metric 

This section investigates the path of the geodesic by considering the direction of its 
tangent vector and the jumping mechanism used in the algorithm. Unless otherwise 
stated, the Einstein summation convention is used on all quantities in component 
form. 


3.1 The Geodesics under Conformal Mapping 

Theorem 3.1 The level curves and the gradient of the objective function </> are 
the attractors of the geodesics on any manifold conformally related to Euclidean 
space. 


3.2 Jumping Mechanism 

Occasionally, the geodesic can be confined in the neighborhood of a local maximum. 
Here we discuss a method to estimate the direction of a neighboring maximum from 
the local maximum using a geodesic so that a jump can be implemented to restart the 
geodesic along that direction. 

Let l be the length of the geodesic, 7. We define the jumping direction to be 
along the vector 

Ax := - J <j)(x)x dx — - j x dx, 

where the integral is over the geodesic and (f> being the normalized <f> over t, . In 
practice, this is approximated by the sum over all steps along the geodesic 

1 T ~ 

^ x - ^ _x *]’ 
t =1 

where T is the total number of steps and x t £ 7. This is just the difference of the 
weighted mean and the mean position vectors along the geodesic. Suppose that a 
neighboring maximum exists and the geodesic is symmetric about a local maximum 
(as it would usually be the case if the geodesic is trapped, for instance, as in Figure 
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[2}. Then the weighted mean would be slightly biased towards the neighboring max¬ 
imum. And so Ax would be pointed towards the neighboring maximum. We used a 
decreasing jump distance for each jump. This is by no means an accurate estimate 
of the direction to the next maximum. However it has been proven to be sufficient 
for our algorithm to discover the global maximum in many objective functions of 
multiple maxima. 


3.3 Solving the Geodesic Equation with a Quadratic Approximation 

In this subsection we discuss the quadratic approximation used to solve the geodesic 
equation iteratively. We give an estimation of the adaptive step sizes to ensure that 
the approximation is valid. The geodesic equation is 

d 2 x i (t) , ni dx j (t) dx k {t) _ n ^ 

^ + 1 > k ^t dT~"- (2) 

In the neighborhood of x t , the (discretized) approximation to the solution of the 
geodesic equation is 

x t+1 = x t + v t 5t + c t (St) 2 , (3) 

where v t is the unit tangent vector of the geodesic at x t , St is the step size, and 

Ct = \~^t = ~ 2 ( Vt ' V< ^( Xt )) Vt ]- (4) 

The tangent vector is estimated by the (normalized) difference x t — x t _i and we 
set the initial tangent vector to be the gradient, v t= i = V<^>(x t= i). Note that the 
quadratic approximation ([3]) is not valid when there exist a component i such that 
0{v t iSt) 0(c t i{St) 2 ), as the approximation is only accurate when the quadratic 
term is small. The value of St when the linear term equals to the quadratic term in 
magnitude is St = tc, where 

tc = min 

i 

If the quadratic term has opposite sign to the linear term, then x l t+1 = x\ when 
St = tc for some component i. Also, when St = \tc, x\+i ~ x t is maximized. 
Therefore, at every time step, the algorithm chooses a step size of St = | tc, if it 
is not smaller than the specified lower bound on t (to be explained in the algorithm 
section). 

Since the geodesic aligns itself with the gradient (or the level curves) as shown 
in the Section [3~i~| in both cases, the step sizes are 

II II - 3 1 

H Xt+1 4 |V0| ‘ 

It can be found by substituting St = ^tc into Equation [ 3 ] and setting Vi = 1 in the 
component parallel to the gradient (or the level curves). Furthermore, as the geodesic 
travels towards a maximum following the gradient, it has a linear rate of convergence 
similar to gradient descent. 


vu 

Cti 
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Fig. 1 Geodesic traversing through multiple local maxima in the search space. 


Fig. 2 Geodesic trapped in a local maximum. 


4 Algorithm 

The algorithm has two parts. The first is a geodesic guided optimization (GEO) al¬ 
gorithm. It estimates the geodesic using the quadratic approximation for a total of T 
steps. The step size is adaptive and bounded by the validity of the quadratic approx¬ 
imation. Quasi-Newton (QN) optimization may be performed, using points along 
the geodesic as the starting points. Figure [I] shows the estimated geodesic moving 
through three local maxima. The algorithm returns the location of the maximum and 
its function value along the geodesic, or the one obtained by QN, whichever is the 
largest. 


4.1 Choosing Parameters 

There are cases where the geodesic fails to reach multiple maxima, for instance, as 
in Figure [2] When the objective function is highly oscillatory, the global maximum is 
less likely to be found by the geodesic. Furthermore, a lower bound on the step size 
St lb must be specified as an input parameter to prevent the step size to become un¬ 
practically small in regions of large gradient, but the choice of an appropriate lower 
bound for any objective function is difficult (if not impossible) to determine. Intu¬ 
itively, a large St lb would allow the geodesic to escape from local fluctuations. On 
the other hand, it may prevent the geodesic from visiting the global maximum. 

The second part of the algorithm. Sequential GEO (SGEO), is developed to 
overcome these difficulties. Information from the geodesic is obtained and passed 
to SGEO. It includes an estimate of the direction of a neighboring local maximum. 
Ax, an indicator, /.;, to denote whether the geodesic is trapped in a local maximum, 
and the average distance between the starting points and the end points of QN, R, to 
determine whether the objective function is oscillatory. 

SGEO calls GEO sequentially with decreasing StLB for N times. In each sub¬ 
sequent run, StLB is reduced by a factor of a , determined by requiring that StLB in 
the last run to be 1000 times smaller than that in the first run. The next GEO run 
starts from a position obtained by translating x* along Ax, with the magnitude and 
method of the jump determined by k. The two mechanisms described above assist 
in the escape from local maxima. In the case of oscillatory functions, QN is not per¬ 
formed, allowing for a larger number of GEO runs. The algorithm first assumes a 
non-oscillatory function, and adaptively adjusts its parameters suitable for an oscil¬ 
latory function by setting a threshold on R. Finally, we impose a stopping criterion 
to improve its computational cost. The technical details are described in the follow 
subsections. The only inputs needed are the number of GEO runs, N, total number 
of steps Nt, the stopping threshold, N t h, which only depends on the dimensionality. 
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and the initial St lb which is only dependent on the dimensionality and the size of 
the search region. 


4.2 The first component, GEO 

GEO estimates the geodesic corresponding to a conformally Euclidean metric with 
the conformal factor given by the objective function up to T steps and evaluates the 
objective function (f>{x t ) at every time step, t, along the geodesic. 

At each t, the normalized tangent vector is used to evaluate x t+1 in Equation 
[5] Since the step size is at most of order ()( L/V,c6) as discussed in Section 
The geodesic tends to get trapped in regions of large gradient. To avoid this problem 
we introduce a lower bound on the step size, St lb, so that the step size is St = 
max{0.5fc) ^lb}- The lower bound also ensures that the geodesic has length of at 
least TSt. 

Now, consider the case where tLB tc ■ The trajectory is dominated by the 
quadratic term c t (5t) 2 . But at t = 1, the tangent vector is the unit gradient and so 
c t= i = — V0/2, the geodesic moves against the gradient. An additional minus sign 
is introduced to c t whenever tus > 0.5 tc at t = 1. A backward geodesic that 
initially moves against the gradient is also estimated by using the initial condition 
v t= i = -V0. 

For both the forward and backward geodesics, Quasi-Newton optimization can 
be performed at every Tqn steps. Whenever x t+ i is outside the search region, the 
algorithm sets x t+ 1 to be a random point sampled uniformly within the search region. 
The following quantities are also evaluated to pass to SGEO, the mean distance be¬ 
tween the Quasi-Newton initial position and the solution R, the normalized mean of 
</>(x t )x t — x t over all t, and an integer k £ {0,1, 2} which characterizes the degree 
of locality of the geodesics, 

! 0 if <Sf is not unique in the forward geodesic within tolerance Vt. 

1 if S)f is unique in only the forward geodesic within tolerance Vt. 

2 if (/>£ is unique in both geodesics within tolerance Vt, 

where S>t is the larger of </>(x t ) and the optimized value with QN. The algorithm 
returns the maximum objective function value and its position along the geodesics as 
well as the information needed to pass onto SGEO. 


3.3 


4.3 Algorithm 1: GEO, Geodesic Guided Optimization 

Matlab’s fminunc() function is used for the Quasi-Newton optimization. Its input 
parameters are set as Maxlter = 200, ToQ = 0.05, Tolx = 0.01, where the last 
two quantities are the tolerances in <j> and x, respectively. 

The set of input parameters for GEO is 

- T, the number of time steps along the geodesic, 

- Tqn, the number of time steps between each QN call. 
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- xo, a vector the initial point of the geodesic, 

- </>(•), the objective function, 

- V(/>(•). the gradient of the objective function, 

- L, a vector containing the lower bounds of the search space, 

- U, a vector containing the upper bounds of the search space, 

- St lb, the lower bound on the step size, 

- sqn, boolean variable denoting whether QN is performed. 

The algorithm is as follows 

1. Forf = 1 : T 

(a) Calculate (f>(x t ) and V</>(x t ). 

(b) If t = 1, set x t := xo, k := 1 and Sx t := V(/>(x t ), 

i. else, set <5x t := x t — x t _i. 

(c) If mod(t, Tqn) = 0 and sqn = 1, 

i. call QN(<f>(-), Xt) and obtain x*}, 

ii. set Rm ■* := | |x* — Xf 11 and set m := m + 1, 

iii. else set <j>*^ := </>(x t ). 

(d) Calculate the normalized tangent vector v t := <5x t /||<5x t ||. 

(e) Calculate c t := |[V0(x t ) - 2(v t • V^(x t ))v t ]. 

ff) Set tc '■= miiij \v t i/c t i\, where * € {1,..., D} denotes the i-th component. 

(g) Set the step size St := ma x{^tc,StLB}- 

(h) If St = St lb and t = 1, change the sign of c t := — c t . 

(i) Calculate x t +i := x t + v t St + c t ( 6 t) 2 . 

(j) If xu + i)j < Li or x/ t+ i)j > Ui for any component i, sample X( t+1 ) from a 
uniform distribution in [L, U]. 

00 Calculate Ax^ := yE[ = i( x t^ (1) - x t)], where cj)* t (1) = 4>* t {1) / £ 4>*t (1) 

(1) Set <^*(0 = max^^l 

(m) If |0*O) — < Tolf(j> *0) for all t, then set k = 1. Else set k = 0. 

2. For the backward geodesic, step (1) is repeated with the following adjustments, 

(a) 5x t is defined as Sx t := —V<j)(x t ), in step (lb), 

(b) := QN(</>(■), x t ) in step (lc(i)), 

(c) Rm := ||x* — x t 11 in step (lc(iii)), and 

(d) omitting step (lh). 

(e) Ax^ := ( x t?t (2) - x t)] in ste P (Ik), 

ff) Set 4 >*( 2 i = max</>(i 2 i in step (11). 

(g) If | (jfk 2 ) — (j)*W\ < Tolf(j)*^ for all t and k = 1, set k = 2. 

3. Set Ax := i[Z\ X W + Z\x( 2 )]. 

4. Return 0* := maxl^ 1 ), c/>*( 2 )}, x* := (x|<^>(x) = R := mean{R(0, R®}, 
Ax := Ax/||Z\x|| and k. 


4.4 The second component, SGEO 

SGEO runs GEO sequentially with different parameters. Let U and L be vectors 
denoting the upper and lower bounds of the search space and let A = min(U — L). 
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This algorithm checks whether the objective function is highly oscillating. We use the 
following criterion that an oscillatory function must satisfy: R < 0.1 A\/D in any the 
first two GEO calls. The reason for limiting to just the first two GEO runs is that St lb 
gets smaller after each consecutive runs and it is more likely for R to be small even 
for non-oscillatory functions. In the case of a high dimensional (D > 10) oscillatory 
function, no Quasi-Newton optimization is performed to allow for a higher number 
of GEO runs. Both of which are crucial in locating the global optimum of highly 
oscillating functions. 

The algorithm uses a procedure similar to annealing to reduce St lb for each 
GEO run. Initially, St BB °^ is set to be AyfD /100, where n denotes the ?r-th GEO 
run. Then the lower bound on St is lowered such that St r l B = a n St B l B °' > ,a € (0,1). 
The factor a is chosen such that St^ B N ^/St^ B °^ = 10 -3 . 

(n) 

After each GEO call, the initial value, Xq , for the next GEO call is estimated 
depending on the value of k passed from GEO. Intuitively, Ax would be a vector 
pointing roughly towards a neighboring maximum. For k = 0, the local geodesic is 
not trapped, 

x<" +1) = x* w + (a ri l)4x ( ’ ,) . 


For k = 1, the forward geodesic is trapped and we set Xq" + 1 ^ to be further away from 


An) 


4" +1) = x* (n) + (a A) A: 


■ (n) 


Finally for k = 2, when both the backward and forward geodesics are trapped, the 
method using Ax becomes ineffective as the objective function has similar values 
along the geodesics - Ax (,, ) points in the same direction as x tn> . Therefore we simply 
set Xq" + 1) to be a point reflected across the midpoint of the search space from Xq H \ 


An+ 1) 


L + U 


M( 


L + U 


A stopping criterion is imposed to reduce the computational cost. Let = 
{(j>*( n ~ 1 ),..., 4 >*(")} be a series of maxima found up to the ??.-th GEO run, <j>* = 
max^( n ) and N* be the number of elements in 0 < "' that are within tolerance of <]/ . 
The algorithm is stopped if N* > N t h(D), where 


f 5 D < 10 

N t h(D) = < 10 10 < U < 20 
[20 20 < .D < 50. 


The algorithm returns rfi* and x* = |x|</>(x) = (j>*}. 


4.5 Algorithm 2: SGEO, Sequential Geodesic Optimization 

Let A = min(U — L), the set of input parameters is 
- N , the number of GEO runs. 
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- Nt, total step number, 

- St^g°\ starting lower bound on the step size, 

- N th (D). 

The algorithm is 

1. Set A := min(U-L), then set {N, a, N T , 5t { ™= 0) ,T QN , s QN } = {20,0.7, 500, Ay/D/100, 10,1}. 

2. Calculate T := and sample x[ | ra_1 ' ) uniformly in [L, U]. 

3. For n = 1 : N 

(a) Calculate StSfy := aSt^g 1 ^. 

(b) Obtain {</.*("), x*("), R, Ax '^, k} by calling GEO(x' n) , 6t ^, s QN ,T, T QN ). 

(c) If D > 10, f? < Q.IAVD, n <2 and sqm = 1, 

i. set {N, a, N T , sqn} := {400,0.98,4000,0} and 

ii. break and restart current loop with the parameters in the above step in 
place of those in step 1. 

(d) If k = 0, set Xq" +1 ' 1 := x*^ + (a n A)Ax \ Else if k = 1, set Xg" +1 ^ := 

x*( n ) + (aA)Ax \ Else set XQ n+1 ^ := + a n A(^A^- — Xq”^). 

(e) If n = 1, set </>* = Else set <f>* := max{</>*, 

(f) Find N*, the number of instances such that \(f>* — <f>*( n < Tolfifi*, n' £ 

{1 ,..., n}. 

(g) If N* >N th (D), exit loop. 

4. Return (f>* and x* := {x|</>(x) = <^>*}. 


5 Numerical Experiments 

In this section we compare SGEO with other algorithms on test functions commonly 
used in the literature. The objective functionsQused have dimensions ranging from 2 
to 50. Sixteen of which are reasonably smooth, the other twelve are oscillatory. All 
calculations are performed in Matlab. 

Table[l]shows the number of failures in locating the global maximum for exist¬ 
ing methods. A success is defined when the estimated function value is within 5% of 
the maximum value. If the maximum value is zero, a success corresponds to finding 
a function value less than 0.05. As it can be seen in Table[l] the Global Search (GS) 
method in Matlab’s global optimization toolbox outperforms Quasi-Newton (QN), 
Trust Region (TR), and Genetic Algorithm (GA). 

In Table [2] we justify the use of QN and the jumping mechanism in SGEO. 
The algorithm without QN fails in high dimensions whereas without jumping the 
algorithm fails in oscillatory cases. 

For finding a global maximum from the chosen test functions, we found that the 
global search method is the best among all commonly used optimization method. We 
therefore did an extensive comparison with the global search method in Table [3] We 

1 These are obtained from http://www.sfu.ca/~ssurjano/optimization.html We have used the log of 
the Hartman and Goldstein-Price functions as the values of these functions vary across five orders of 
magnitude within the search space. 
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found that SGEO can discover the true global maximum with a higher chance than 
GS in the objective functions tested. Our method is more accurate than the global 
search in many test functions. There is only one function that our method is not as 
good as the global search. At the same time, the computational cost represented by the 
number of function calls and the computational time remains similar in most cases. 


Table 1 Number of failures over 100 runs for Quasi-Newton (QN), Trust Region (TR), Genetic Algorithm 
(GA), and Global Search (GS). 


Table 2 Number of failures over 50 runs for SGEO, without Quasi-Newton (with T = 200), and without 
jumping. 


Table 3 The number of failures (Nf a u ure ), computational time, and the number of function calls (N ca u) 
over 50 runs of SGEO and GS. 


6 Conclusion 

A new algorithm is proposed in order to find multiple optima of a continuous objec¬ 
tive function. The path constructed by the algorithm follows either a geodesic or a 
contour line. Conformal mapping and the Newton Raphston algorithm are employed 
to enhance computational efficiency. A built-in jumping mechanism also directs the 
proposed algorithm to a more promising search area. A stopping criterion is imple¬ 
mented if the same maximum is found too many times. We are extending this algo¬ 
rithm to handle optimization in high dimensions with contestation. 
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